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Introduction 


In recognition of the advanced design techniques incorporated in the logic ot 
the 6600, the development of the 6600 programming systems included a maflie- 
matical investigation of the library functions intended for the 6600. Previous 
methods of approximations were analyzed in terms of the 6600 and, where 
appropriate, new algorithms were programmed and tested asing a CONTROL 
DATA 1604A in double-precision. 

As a general goal, approximations of the required accuracy were sought which 
were as short as possible with the constraint that the range reduction mechanism 
be of reasonable complexity. In most cases, continued fraction forms appear best, 
although polynomial and rational forms were deri%'ed for many of die functions 
in order not to prejudge during the development phase the form most efficient in 
practice. Therefore, rather than limiting consideration to algorithms optimized 
for the 6600, this document includes a general discussion of methods along with 
presentation of a variety of algorithms selected to meet the general requirements. 
Consequendy, the algorithms and die given coefficients, which were taken di¬ 
rectly from die 1604A double precision results, are designed to achieve the stipu¬ 
lated accuracy but do not take into account the 6600 round-off effects. 

Although the information derived here was gathered during the study of 6600 
library function algorithms, the exact techniques implemented in the original 
library version of the 6600 systems or subsequent revisions thereto are not neces¬ 
sarily contained in this document. 
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Characteristics of the SSOO* 


The CONTROL DATA 6600 is a large-scale, solid- 
state, general-purpose digital computing system 
composed of eleven independent computers. Ten of 
these are peripheral and control processors, each 
with a 12-bit 4096-word memory. The eleventh com¬ 
puter, the central processor, is a high-speed arithme¬ 
tic device. The common element between these 
computers is a random-access central memory of 
131,072 words (of length 60-bits) organized in 32 
banks of 4096 words each. It is the central processor 
whose characteristics are to be considered in select¬ 
ing appropriate algorithms. 

The central processor has ten independent arith¬ 
metic and logical units which operate concurrently 
in the solution of a problem. In addition, it has 24 
operating registers for functional units and 8 tran¬ 
sistor registers for servicing functional units. A 
word length is 60 bits, 48 bits of which determine 
the integer coefficient, 11 bits the biased exponent, 
and 1 bit, the coefficient sign. Execution times for 
floating-point operations are as follows: 

Floating-point add: 

4 minor cycles «4{1(K))(10'*)SECS=4(10‘^}SECS 
Floatinfi-point multiply: 

10 minor cycles = 10(100)(10‘")SECS = 10-" SECS 
•(See Ref. 17) 


Floating-point divide: 

29 minor cycles=29(iOO}(10‘®)SECS 
«2.9(10")SECS 

Hence, relative speeds are 
M=2.5A and D=2.9M*7-25A 

Central processor instructions are sent automatic* 
ally and in the original sequence to the instruction 
slack which holds up to 32 instructions. A branch 
to another area of the program voids the old in¬ 
structions in the registers and brings in new ones. 
Branch orders of the type “GO TO K IF Bi<Bj 
require 6(10*^) SECS and an additional 5(10-0 
SECS for a branch to an instruction which is out 
of the stack. High speed in the central processor 
depends upon minimizing memory references and 
waiting time for unrelated instructions and partial 
answers. 

Since arithmetic computations are extremely effi¬ 
cient and several operations may be done simulta¬ 
neously, algorithms which break down into inde¬ 
pendent blocks which can be computed in parallel 
are desirable. On the other hand, it is surmised that 
division of the interval of definition of a variable 
into many small sub-intervals, requiring the com¬ 
puter to do much testing and branching to other 
blocks of the program not in the stacker, is not very 
efficient for the 6600. 
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Method of Testing 


Let F(x) be the function to be considered as the 
correct one, and let Y{x) be the approximation to 
F(x) being tested. Define A(x)=| F(x)—Y(x) | as 
the absolute error in Y over some x range, and 


as the relative error in Y over this 


I F(x) I 

range. For floating-point subroutines, accuracy is 
defined by the number of first correct significant 
digits, so that if 


the n first significant digits arc correct. For fixed- 
point subroutines, the absolute error measures 
accuracy. If 

A < S { 10-^“**0 


then n digits after the decimal point are correct. 
Since all approximations considered here are to be 
in floating-point, it shall be required that 

R<2'^^-'L775(10-^^), 

which is the basic roundoff due to the size of the 
6600 register. 

The algorithms Y(x) described here have been pro¬ 
grammed on the 1604 in FORTRAN 63 using 
double-precision, and compared with double-pre¬ 
cision library routines F(x) supplied by Palo Alto. 
This provided a test for truncation error in the 
algorithm itself, but, of course, gave no effect of 
6G00 roundoff since the floating-point word length 
for the 6600 in single-precision is between those of 
the 1604 used in single and double precision modes, 
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III. General Methods of Approximation* 


Three forms of approximations have been consid¬ 
ered—(1) polynomials, (2) rationaIs» and (3) (trun¬ 
cated) continued fractions, together with techniques 
for improving convergence and for converting from 
one form to another. 

The following sections attempt to give in capsule 
form the theoretical background for the approxima¬ 
tions developed in later sections, though not all 
methods discussed were actually used in these 
approximations. 


T„.i(x)=2xTo(x)- T>-,(x) 
with Ti,(x) = 1 and Tj(x) = x. 

Writing 

T (x)=C« +C'x+C2x^ + ... + C“x" 

n' ' 6 n D n 

the coefficients are computed from 
C~ = 0 if (ti'^m) is odd 


A. POLYNOMIAL APPROXIMATIONS 

1 . Truncated Taylor-Maclaurin Power 
Series 

One of the most common approximations upon 
which many others are based is the Taylor- 
Maclaurin power series truncated to m terms, 

f(x)^2-, I X I <a, where m is 

Dsu nl 

chosen sufficiently large to insure the desired accu¬ 
racy. The absolute error is less than the maximum 
of the absolute value of the first neglected term 
over |x|<a. Such approximations usually require 
too many terms to be used directly unless the inter¬ 
val ; X {<a is so subdivided that the full number of 
terms is used only a small proportion of the time. 
(In this connection, see Ref. 13 on the subject of 
‘'Partitioned Polynomials’). More often, Taylor- 
Maclaurin series provide a starting point from 
which more efficient routines can be built. 


2. Chebyshev Expansions (See Ref. 3 and 5) 

Denote by Tn(x)scos n9 the Chebyshec' polyno¬ 
mial of degree n in x=cos 0. It is simply the poly¬ 
nomial in X obtained by expressing cos in terms 
of cos 0, then replacing cos $ by x. Note that 
: X I <1 and that | Tn(x) <1. The polynomial Ta.i(x) 
attains its greatest abso ute value, one, in the inter¬ 
val [ — 1,1] at n + 2 points (including endpoints) 
with alternating sign. These polynomials may be 
generated by the recursion 

•(See Ref, 5, 4. 1,2 and 18) 


2 (-i) ‘ 

r/(n + mj is even 

Alternately, the shifted Chebyshev polynomials 
T;(x) = T„(2x-1), 


1-hcos 6 


i.e., 


obtained by the transformation x = 

by replacing x by (2x—1) in the polynomial Tn(x), 
can be used for the range 0<x<l. Coefficients for 
T^(x) arc found from 

c-2~ [2(::) - (“r)] (-1)"“' 


Tables of Chebyshev coefficients may be found in 
Ref.s. 3 and 13, and also in publications by the 
National Bureau of Standards. Reference 7 pro¬ 
vides an excellent source of information for locat¬ 
ing appropriate tables. 

The Chebyshev expansion of f(ax) over the inter¬ 
val — a<ax<a, truncated after the mth term, is 
written 


f(ax)^ +V C„(a)T.(x),]x|<l 


where 




C.(a)= ^ / f(ax) T„(x)(l-x^) dx. 
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The tnincation error is approximately C,B^,(a) 
Tm*i(x). Since the coefficients in this expansion 
depend upon values of f{ax) in the entire interval 
( — a, a) rather than only upon values at zero, die 
approximation is more efficient and converges 
more rapidly than the Taylor^Maclaurin series as 
n increases. However, the coefficients Co(a) are 
often extremely troublesome to compute accurately. 

3. Telescoped Polynomials (See Ref. 3 and 5) 

If f(x) is an arbitrary polynomial of degree n+1, 
the "best” polynomial approximation of degree n 
in [-1,1] is p8(x)*f(x)^a„4jTn*i(x) where a^*! is 
a constant chosen so that the coefficient of x*** on 
the right-hand side vanishes. In the sense of the 
next Section (III.A.4,) the polynomial Tn.,(x)/2** is 
the unique "best” approximation to zero in [—1,1] 
with leading term exactly 


or, X 




To*x(x) 


h the unique "best” approximation to x“** of de^ 
gree n. This fact forms the basis for a telescoping 
procedure described in the following manner. 

Let f(x)aao+a,x+a8X* + ... + a.*iX"*\ |xi<l or 
xc [0,1], represent f(x) to the required accuracy in 
the interval Usually f(x) is represented by a trun¬ 
cated Taylor-Maclaurin power series. If Ae range 
of X is not [-1,1] or [0,1], a simple transforma¬ 
tion can be made to a variable with such a range. 


Now 


or 


Cn*i 






where 


t “C"* /C"*’ 

ni B-'l 


Next try replacing x"** by 

T^Jx) 

, i.e., substitute 

n+1 


T (x) 

t+l' ' 




B+l 


“t„X“-ta-iX"-*- . . . -tiX-t«. 


Letting 

a'=a^ — a^^^ t. for i = 0, n, the result is 


f{x)s a;; + a^xH-+... + a'x" + 


Since 


Cd*i 

n+1 


Tn4i(x) |<1, wehave 


a T ,(x) 

B«1 a+1' ' 

< 







B + l 


= E,. 


Let Ef, be the original truncation error incurred by 
iLsing terms of the Taylor series only through 
(n + 1), and let E be the allowable error. 


Then if Eo+Ei<E, the term an-i To*i(x)/Coti 

may be dropped, and the result is a polynomial 

approximation for f(x) of degree n with error less 

than E. The process may be repeated so long as 

S Ei<E. When using Tb4i(x) for telescoping, the 
1*0 

highest power of x will decrease by two with each 
substitution, since To«i(x) contains only alternate 
powers of X. howes'er, contains all powers 

of X from 0 through n + 1, so that its use in tele¬ 
scoping reduces the highest power of x by one with 
each substitution. Note that telescoping a power 
series using Chebyshev polynomials is not equiva¬ 
lent to the Chebyshev expan.sion described in Sec¬ 
tion 1ILA.2. (See Ref. 5, p. 12 for an example). 

4. Best-Fit Polynomials (See Ref. 22 and 4) 

Let f(x) be a given function continuous on the 
interval [a, b], and let g(x) be a given weight 
function, continuous and positive on [a,b]. That 
polynomial, P^(x) for which 


Max 

(». b] 


p;(x)-f(x) 

g(x) 


— Min 


is called the Chebyshev-approximant or best-fit 
pol>Tiomial of degree n (in the sense of Chebyshev) 
with respect to the weight function g(x). The 
weight function allows the option of minimizing 
either absolute or relative error by taking g(x)^l 
or g(x)=| f(x) I respectively. 
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Let f(x) be an arbitrar>' single-valued function de¬ 
fined in the closed interval [a, b] and let Pq(x) be 
a polynomial of degree n such that the deviation 

fo(x)=f(x)-p,{x) 

attains its greatest absolute value L at not less than 
n+2 distinct points in fa, b] and is alternately +L 
and — L at the successive points. Then pn(x) is the 
best polynomial approximation of degree n to f(x) 
in fa^b] in the sense that the maximum value of 
I rn(x) I is as small as possible. 

These conditions imply a set of (2n+2) equations 
for L, the (n + 1) coefficients of p,(x), and the n 
critical points, Xj, at which the value ±L is at¬ 
tained (the other two are endpoints). The assump¬ 
tion that «'^(x) is continuous on [a. b] and is zero 
at each x» is usually required. Hence the 2n+2 
equations to be solved are 

«D(Xi)=f(Xj)-pn(X|)»{-l)*L 


is known for the best polynomial approximation of 
given degree to an arbitrary function f(x). 

B. RATIONAL APPROXIMATIONS- 
QUOTIENT OF TWO POLYNOMIALS 

1. Pade Approximants and Table 

(See Refs. 6,5) 

eo 

Given a power series P(x) = 2 and a pair of 

1=0 

non-negative integers (m, n), there exists a uniquely 
determined rational fraction, R"(x), whose numer¬ 
ator and denominator are of degrees less than or 
equal to m and n respectively and whose expansion 
in ascending powers of x agrees term by term with 
P(x) for more terms than that of any other such 
rational fraction. R“{x) is called a Pad6 approxi- 
mant of P(x) and the associated table formed by 
putting R*(x) in the (n + l)st row and {m + l)st, 
column, n, m=0, 1, 2,..., is called a Pad^ table. 

Let 


for 1 = 0 , n+1 




for i= 1, n . 


P(x) = 2 Cix'. A„(x) = 2 aix'. 


B„(x) = 2b.x-. 

i*9 


Their solution requires some sort of iterative pro¬ 
cedure. (See Ref. 22). 

If Pb(x) exists, it Is unique. If f(x) is a c<mtinuous 
function in [a, b], then there exists a unique poly¬ 
nomial of best approximation of given degree. If 
f(x) is a polynomial of degree n + 1, the best poly¬ 
nomial approximation of degree n in [-1,1] is 
Pii{x) = f{x)-an.j T„,i(x), where an*i is a constant 
chosen so that the coefficient of x"*‘ on the right- 
hand side vanishes. No simple explicit expression 


Then there are m+n+2 coefficients ai, bi to be 
found. 

Hence, if we write 


Bl(x) 


2 Cux** 

krii«ai*i 


P(x)B„(x)-A„{x)*B„(x) 2 C.x'-s 2 d,x^ 


C„bo + (C,bo+C„b,)x + ... + 2 c..,b.x>' + ... + 2 Co-,b,x" + 


2 C..,-,b,x“*' + ... + 2 Q-.b.x" + ... + 2 + 

i*0 


+ ... - fao+a.x + ... + a„x"] = + — 
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and set all coefficients of x through x“*“ to zero, 
the result is 

2 Ck-ibi —a^—0 for k^O, n 

UU 

n 

2 Ck-ib| —aii=0 for k=n + l, n+m 

I 70 

ak = 0 for k>m 

which provide n+m+1 equations to be solved for 
a I and bi. These equations will not contain ais for 
k®m+l, m+n, so will provide n equations for 
finding (n + 1) b|S. Since one of the bi’s must be 
arbitrary, choose bn® 1. 

Thus, to solve the above equations for the coeffi* 
cients of the rational approximation 

set 

bc = l 

and solve the n equations, 

n 

2 Ck-ibi =0 for k—m+1, m+n for bi, i=l, n. 

I 0 

Then solve the (m4' 1) equations 

miB(k.i)) 

2 Ck-(bi= ak, k=s0. m 

Isg 

for ai, i*0, m 
and 

Pad6 approximants are most useful for m®n or 
m=n + l. P(x) is normally taken to be the Taylor- 
\laclaurin expansion of the function to be approxi¬ 
mated and n and m chosen so that the terms of the 
series through yield an approximation with the 


accuracy desired. If m=n, the resulting rational 
function is more accurate than the series through 
x“*®; i.e., tlie rational function actually agrees with 
more terms of the series than required. An error 
estimate is gi\'en in Ref. 5, p. 14, which requires the 
evaluation of the quotient of two determinants of 
orders (n+l)x(n+l) and nxn for the case m=n. 
The elements of the determinants are Ci s. 

2. Maehly's Method (See Rets. 5 and 13) 

Rational approximations may also be derived from 
Chebyshev expansions as shown by H. Maehly. If 
a function f(ax) in the inten^al — a<ax<a is ex¬ 
panded in its Chebyshev' series, 

f(ax) = 2 CiT|(x), (See Section III.A.2.) 

and 

ni n 

A„(x)=2aJ*(x). B,(x)»2bT,(x) 

• Ttt . I 

then the (n + m+1) unknowns ai, bi (bo®l) are 
determined from 

(2biT,(x)) (2C,T,(x)) - (2a,T,{x)) 

= S h.T.(x). 

Using the relation T™*„(x)+Tm-n(x) = 2Tm(x)Tn{x), 
the following system of n+m+1 linear equations 
is obtained for the z, and ht: 

ao=Co+ — 2b|Ci 
2 1.1 

a.=C, + -^+l2:b. (C,.,+C ,_,.),] = l.m+n 

where aj=0 for ]>m and bj=0 for j>n. Note that 
it is nec*essary to find the coefficients Ci before 
Maehly's Method may be used. An error estimate 
is given in Ref. 5, p. 16, and a detailed example 
from Arcs in x is given in Ref. 13, pp. 123-131. 
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3. Best-Fit Rational Approximations 

{See Ref. 22) 

Let 

m 

Aen(x) = 2 be the numerator and 

l-u 

n 

Bn(x) = 2 biX' the denominator of a rational 

1 

function RN(x)*An,(x)/Bn(x), N = n+m. 

In the .same manner as was done for polynomials, 
a rational best-fit approximation of order N to the 
continuous function f(x) on the interval [a, b] is 
defined as that rational function R* (x) for which 

max ----= min. 

g(x) 

where g(x) is a given weight function, continuous 
and positive in [a, b], e.g., g{x)*l or g(x)s| f(x) |. 

N, Achieser has shown that R* (x) is uniquely char- 
acterijted by its error curve, 


R;(x)-f(x) 



assuming its maximum absolute value sufficiently 
often with alternating signs. Arguments x* for 
which the maximum absolute value is assumed are 
called critical points, 

An error curv’e has standard form if it meets the 
additional requirements that it has exactly N+2 
critical points, the first and last of which are end¬ 
points of the interval, and has a continuous deriva¬ 
tive with respect to x which vanishes at the critical 
points. 

If an error curve has standard form, it is necessarily 
the optimal error curv’e corresponding to the best- 
fit rational R* (x). The converse is not necessarily 
tnie—the optimal error curve need not have stand¬ 
ard form. 


8*(x) in standard form yields 2X + 2 equations in 
the 2N+2 unknowns x’, i=l, N; ai, i=0, m; 
bj, i®l, n, and the maximum error, Since one of 
the coefficients in R^.(x) is arbitrary, assume b(,=l. 

8(x0=(-l)‘X,i=0.N+l 

a'(x,)=o,i= i,N 

This is a non-linear system of equations whose 
solution re<iuires an iterative procedure. The term 
“direct’* method is used to indicate that the coeffi¬ 
cients of the hest-fit rational are computed directly 
from the equations above, whereas an “indirect 
method determines the corrections necessary to 
modify a fixed approximant (e.g., a Fade approxi- 
mant) to obtain the best-fit rational. Details for 
several direct, indirect and combined methods, due 
largely to the late H. Maehly, may be found in 
Reference 22. 


4. Conversion of Rational Approximations 
to Continued Fractions 

Any rational fraction may be converted to an 
equivalent wntinned fraction. The continued frac¬ 
tion form evaluated from bottom to top (See Sec¬ 
tion III.C.I.) is nearly always more efficient (i.e,, 
requires fewer operations) to evaluate than the 
rational form. 

Assuming the degree m of the numerator is ^ the 
degree n of the denominator, set the rational func¬ 
tion identically equal to the continued fraction: 

a>i+a|2-fa^2^ -f- ... 4- amZ^ 
bo+biZ-l-bsZ* + ... + boZ" 

sC+ ^ 

" (z+B0 + (z+B 2 )+ ... (z + B„)* 

Since the a, and bj are known, the Ci and Bi are 
found by equating coefficients of like powers of z 
after converting the right-hand side to a rational 
form. The resulting equations are non-linear, but 
easy to solve. Results for n=2,3,4 are given in Ap¬ 
pendix A as COF' 2, COF 3 and COF 4 respectively. 
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C. CONTINUED FRACTIONS 

(See References 6 and 8) 

1. Notation and Methods of Evaluating 


2. Error Estimate 

In determining the size of n required to give a 
desired accuracy, the following “determinant for¬ 
mula” is often helpful: 


Let 


F = bo+ai 

bi+a* 

bj+ai 

b.i+ ... 


be a continued fraction. 

F may also be written in either of the following 
two notations: 


F=bo + 


^ 

bi+ b3+ ... 


An _ An-i > ^n 

B. B.-i Bo-, Be 

where Bp., and B,, are defined in die previous 
Section 1. Note that this is merely an expression for 
the difference between the (n—l)st and the nth 
approximants and says nothing about F—Fb. No 
simple formula for the error F-Fa is known except 
for special cases. 


3. Equivalent Continued Fractions 

The continued fraction 


F=bo + 





F=C„bo+ 


CuC^a, C}C:a3 

C,bi + Csb- + 


Crt-,Cnan 
Cnbn + 


Assume in the serjiiel that no bi, i>l is zero and where C, # O.has approximant En/Dn where 

that the continued fraction converges. F can be 

evaluated for n terms in several wavs: En“C(,Ci.. .CnA, 


(a) Top to bottom 


An ba An-j+aftAn-j | 

Bn=b„ I 


n = l,2,... 


D„ —C, . . . CfiBn 
Dn® Bo® 1 

and An/Ba is the nth approximant of F as defined 
in 1. 


where A , = 1, Ao = bf.B i=0, Bg = l. 


Ff,® An/Bp is the nth approximant to F. 


We have 

D„ B„ 

\ 

and both coni erge or diverge together. 


(fe) Bottom to top 


Pt=aD 


Qi®bB 

Pi Qi -1 I 

I i=2,n 

Qi®bn-t-i Ql-i + Pj.j 1 

F„ = b„+ ^ 

Q. 


If Co=l, the two continued fractions F and F are 
equivalent. This equivalence is useful for chang¬ 
ing the form of a given continued fraction, say to 
one in which the numerators (or denominators) are 
all one. 


4. Even and Odd Part Contractions 


The even part of a continued fraction is the con¬ 


tinued fraction whose approximants are 



simi- 
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larly, the odd part 15 the continued fraction whose 5. Functions Expressed ds Continued 
approximants are Fractions 


^2n*\ 


.n=0.1.... 


The even and odd part contractions converge if the 
original fraction does and to the same value. 


Each of the previous Sections 1 through 5 applies 
when the at or b, are either constant or a function 
of X. It is assumed that the continued fraction F{x) 
converges for j x j<€. 


Fev«n —bo + 


ajbs 


(bjbs + aj)" 
a4asb;bd 


(b2ba+a3)b4+b2a« 

aea?b4b6 


(b4b5+a5)be + b4a« — (b8b7 + a;)b»+b6as“ 


b^bx + a j aidsb^/b i 

bi {bjbg + ajjbs+bjBs — 

a^a^bibs a^a^b^b; 

(b8b4+a4)b5+bsa5 — {b5b4+a«)b7+bsa7 — ... 


Especially useful are the continued fraction expan¬ 
sions of Gauss for the functions tan x> arctan x> e’ 

and loge 73 “ • These are reproduced in Appen- 

dix B together with some of their successive 
approximants» Ar/B^. 


As an illustration of the uses of some of the rela¬ 
tions given in this Section C, consider the Gatissian 
continued fraction for e* (Appendix B,3.a.): 


5« Continued Fractions Equivalent 
to Series 


X X X X X X X 

1- 2+ 3- 2+ 5-2+7- 


The series C.i + Ci+ ... +Cn... and the continued 
fraction 


Ct c,/c, 

1 - ( 1 +C 2 /C,) - 

C,/C, 

(1 + C,/Cd- {1+C./C,.t) 


The even part c'ontraction of this continued frac¬ 
tion (III.C.4.) is 

, , 2x 2x* 

** (2-x)+ [(6+x)2-2x] + 


are equivalent in the sense that 



4x= 4x» 

[(10+x)2-2x] + [(14+x)2-2x] +••• 


Similarly, the power series C„+C,x+C 2 X=+... 
and the continued fraction. 


2x 2x’ 4x= 4x‘' 

(2-x) + 12 T 20T 28+ 


^ C,x (Cj/COx 

1- 1+(Q/Ci)x - ... 

are equivalent. 

Another method for obtaining a continued-fraction 
expansion from a function defined by a power 
series is the Quotient-Difference algorithm of 
Rutishauser described in References 18 and 19. 


which IS equivalent by Section 111.0.3. to 

2x X* x* x2 
(2-x) + ~6T loT 14 + 

This is now form 3.b. of Appendix B, 
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7. Telescoping Procedures for Continued 
Fractions (Maehly) (See Reference 21) 

In III.A-3. a method was described for telescoping 
a truncated power series by use of Chebyshev poly¬ 
nomials. Maehly has derived a method by which 
the (n+l)st approximant of a continued fraction 
may be telescoped one step to a corrected nth 
approximant. 


To.i(u) is the Chebyshev polynomial of degree 
(n + 1). 

For the uncorrected nth approximant, Rr(x), it is 
known that 


Rn(eu)-f(eu) 




Let 


where 


f(x)= I 

bo bi bs 


+ 






bo i:o bibi»i 


be a convergent continued fraction representation 
of f(x) whose (n+l)st approximant approximates 
f(x) to within the desired accuracy in the interval 
I * I(n+ l)st approximant is 

D / \ I I I l_ I (x) 

Rn.>(x)«= r-' + -j—' + ... + -— = - , ■ 

I bo I b, I b,M B«-,(x) 

where An.i and Bn*j are defined recursively as 
A|,. 1 “ bn. I An ^ An-i I 

/ n>l 

Bn-1 “bti.iBn+an-iX Bn-i ? 

A<i"<iij, At^ofribi 


The corresponding limit relation for R* is: 

R:(eu)-f{eu) 
lim- 

Hence the corrected approximant R^{x) yields an 
“almost best-fir rational. 

Since many of the functions of interest arc "even” 
or "odd" functions, definitions and telescoping for¬ 
mulae will be given for these special cases. First 
replace x by x* in the recursive definitions for Aa^j 
and Ba.,. 


B,»=bn, Bi®b„bi+atX 


Now we may alter either the an, or the bj<, or the 
An and Bn. Formulae are given in Reference 21 for 
all of these cases, but obly those for altering Aa 
and Br, to obtain a new (called R*) are given 
here. 


R* * 



B* 

n 


Aa + y, + X S 

k=2 


Bn + Xyt+X S yi, Bfc-2 

k«2 


where 


Even Fiincriorw; 

An even function, g(x), is one of the form 

g(x)= + ... + + ... for I x 1 <€. 

|b. lb. |b. I I- 

Let 

X = fu where | ti | < 1 and 

Si="'2i(u)= —57^^ = 2 S!i,u’\ whereT 2 „*j(u) 

* k.o 




yk= n ,xs=fU,|u|<l, 


is the Chebyshev polynomial of degree 2n+2. 


and 


Define 




Ta.,(u) 

2 “ 


0*1 




10 




Then 


Then 



Aft+yo+3C" 2 YkAii^i 

ks2 


Bn + X'yi + X- 2 ykBji-2 
k=2 



Ao + yo + X^ 2 ykAn-2 

k=2 


Bp+x^yi+x^ 1 yifBh-2 

ks 2 


and 


lim R;M-g(£u) 


= C„., S<=-‘«> (u). 


Odd Functions: 

An odd function, f(x), is one of the form 


^ bn bx 


lb« 


Lot X—pu where I u 1 <1 and 


T,n-3(u) _ 

k:l> 




25fl*2 

where T„, 3 {u) is the Chebyshev polynomial of de¬ 
gree 2n + 3. 

Define 


yk= -I Sjk.i 1 Tic fork=0.n. 

t k t)| 


and 


lim R;(ru)-f(fiu) 


D. RANGE REDUCTION 

Use of a single approximation to a function f(x) 
over its entire range of definition is usually not 
feasible. Therefore, it is usual to subdivide the x 
range into intervals small enough to provide the 
desired accuracy with algorithms of reasonable 
length, but not into so many subintervals that the 
mechanism for deciding into which interval x falls 
and for acting accwdingly becomes cumbersome 
and time consuming. Depending upon the charac¬ 
teristics of the computer, some balance must be 
struck between the number of intervals of subdivi¬ 
sion and the number of terms (operations) in the 
algorithm. 
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IV. Approximations Obtained for 

Library Functions* 


The library functions considered are square root, 
cube root, sin u, tan u, arctan u, arcsin u, and 
log^u. Other functions may be computed in terms 
of these* Numerical coefficients are not given for 
all of the approximations tested, but those not 
given may be found in SSD Memos, Refs. 9 
through 12. 


A, SQUARE ROOT 

(See References 5,14,15 and 16) 

1. Reduction of Range 

To find VN» N>0, first reduce N to the form 
’x where j <x<l (or j <x<l) and m is 

zero or a positive or negative integer. If this repre¬ 
sentation is to be unique, only one of the two end¬ 
points 1 should be included in the range of x. 
For example 16=2^ • 1=2" ( . 

Hence ^2“' • \/\. 

2. Computation of \/x for 

\ <x<lori<x<l 
4 4 

Vx is computed via a Newton-Raphson iteration 
starting with a first guess, y,, for \/x. Successive 
approximations are found from 

yM = ^(yi +^^ 1 = 1 , 2 ,... 

iterating until |yi.i—Xi |<2'*®. The number of 
iterations required will depend upon the accuracy 
of the first guess, y,. For the various estimates for 
y’l which follow, maximum and minimum number 
of iterations are given for values of N ranging from 
.1 to 10 at intervals of .1. 


•See Ref^ences 1, 2, 5,13,14. 15, 16 and 20. 
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3, Estimates for y. 

a. Best-fit rational derived by Maehly (Ref. 14) 
for y, with max. relative error <2.6(10'*). 

^ b 

Vi = a + —— 
c+x 

a = 3,090315520/>/2 
b * -8.550050013/2\/2 
c = 3.090315520/2 

Max. number of iterations = 4 
Min. number of iterations * 2 

b. Fade approximation in continued fraction 
form (Ref. 5) with max. relative error <2.3(10-*) 
for yi. 


25 

r 5000 / , 15\ 

343 (^"^49/ 


1 

II 



400 


LV 49/ 

V 49^ 

2401 


Max. number of iterations » 3 
Min. number of iterations ® 2 

c. Fade approximation in continued fraction 
form with split range (Ref. 5) and max. relative 
error <10"'in y^. 

For 7 <x< i 
4-2 


Yi ® 


5\/70 

14 


50V70 

49 


(-n) 


(’‘+r4)('‘'"5)‘r9 


For - <x<l 


yi = 


5\/35 


200\/35 

49 


(- 1 ) 




Max. number of iterations = 3 
Min. number of iterations = 1 



The approximation c. gave no better results than 
b. and has the additional disadvantage of using a 
split range. Method a. needs one more iteration 
than b. so that a. takes the time-equivalent of 3 
itera tions + 2D +1M+3 A = 3 iterations + 20A, 
while b. takes 3 iterations +1D+2M+4A—3 iter¬ 
ations + 16.25A, assuming all fractions in a. and 
b, are precomputed. Hence b. is slightly faster, but 
requires 3 more constants than a. Either a. or b. is 
preferable to c. 

B. CUBE ROOT (References 16 and 20d.) 

1. Reduction of Range 

To find vN for N>0> write where 

-<x<l and where k and n are either zero or 

A 

integers with the same sign; k is restricted to the 
values k*0, dil, ^2. Then 


Maximum number of iterations » 5 

Minimum number of iterations = 5 

b. Rational approximation for y» with absolute 
error of about 9(10 at x=l (Reference 20.d.). 

ai 

yi = a.- 

x-f bi 

ao = 1.78781 

a, = 1.91548 

b, = 1.42856 

Maximum number of iterations = 4 

Miniiniiin number of iterations = 3 

c, Continued fraction approximations for yj with 
absolute error of about 9(10 '*) at x=l (Refer¬ 
ence 20.d.). 


2, Computation of v'xfor | ^x<l. 

Find v^x by means of a Newton-Raphson iteration 

starting with a first approximation, yt, for \/x. 
Successive approxi mants are found from 

X 

+ — 

2y: 

iterating until 




»a — a2 

(x+b,)- (x+bs) 

- _ ai(x+b2) 

(x+bi) (x+b2)-a2 

a» = 2.502926 
ai * 8.045125 
h, s 4.612244 
a, = .3598496 
« .3877552 


y,.t-y. I <2‘*^ 

The maximum and minimum number of iterations 
for arguments ranging from .1 to 10 at increments 
of ,l are given for each of the estimates for yi 
which follow. 

3. Estimates for y, 

a. Linear approximation for y, has absolute error 
of about 8{10'*) for x = l (Reference 16). 

yi = A + Bx 
A = .5914(62048 
B = .3319149488 


Maximum number of iterations = 3 
Minimum number of iterations = 3 

Method a. takes 2 iterations more than c., while b. 
takes only one more. Hence, in addition to the 
three iterations required by aU three methods, the 
time-equivalent of 

7M + 2D + 3A=35A is needed for a., 

3M+2D + 3A=25A forb.. and 
2M + ID+4A = 16.25A for c. 

Thus, even if no operations could be done in 
parallel, c. is more efficient than b., which is more 
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efficient than a. Of course, c. would be even more 
effective if parallel operations are done. 

C. SIN u 

1. Range Reduction 

Assume first that the argument, u, has been re* 

duced to the range — ^ T* Th^^ 

2 2 

ranges are considered, 


x|<-and|x[<-g 


SO that for 


< 7T 

— 

- 2 


For IXII < — , let x®u and sin u=sm x 

" 2 


h. 1 X1 < - 
I I - g 


For I u I < — Jet x-ii and sin u*sin x 
• ” 6 


(ir 


n=10, -2.6(10-“)<2- 

i 


and for 




n-7, — _10-’«'^’<2-« 

15! 


3. Telescoped Polynomials (SeeIli.A.3.) 

I I ^ ^ 


Sin 2 ^ ^ with error<2.6{10“'*). 

(2n+l)! 


Lftx»|y .|yl<l- 


For i n I > — , 

0 

let X * and sinu * sin x (3-4 sin* x). 
«3 


2. Taylor-Maclaurin Series {Seelll.A.i.) 

o- _ V 

(2n + l)l 


p^j. I X 1 < — the terms through n=9 provide an 

approximation with absolute and relative errors 

<2-«, while for 1 X [ < -^, n = 6 is sufficient. These 

two polynomials form the basis for much better 
approximations to be developed in the next sections. 


(iy-' 


R.E. < - < -r for I X I < -r, 

“ {2n+l)lsinx (2n+l)! ' '"2 


(r 


R.E. < -- for 1 X I < - 

“ (2n+l)l 1 '6 


-E 


Z-I {2n+l)l ^ 

1180 ' 

--..■(r (i)'" 


(2n+l)! 


Substitute for y'» in terms of the Chebyshev poly¬ 
nomial Ti,(y), 


’19 19 

— y‘‘ - — y*’’ + ... + 
4 y 2 

— y + 
2'* 

^ r'MPI yid*l J- 



~ 1.66(10-'''), 

19! 2‘® 19!2'» 
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and this error added to the original error of 
2.6(10-*") is still <2•^^ so that the term in Ti8(y) 
may be dropped. The result is 


sm X 


(ij 


(2n+l)! 


. C^iioy 

19! 


where is the coefficient of in tl»e 

Chebysliev polynomial of degree 19. Similarly, 
after suhstitiition for y'' in terms of T,:(y), it is 
found that the coefficient of Ti:(y) is about 9(10'*^) 
with the total error still less than 2'^''. so that the 
term containing Tj:(y) can be dropped. A further 
substitution for y'*' results in a coefficient for Tu(y) 
of about 7(10''’’) which is too large to be dropped. 
Tlms we must be satisfied with substitutions for y*" 
and y‘^ only, reducing the degree of the polynomial 
in y from 19 to 15. 

NW, transforming back to the variable x by setting 
y=2x/7r, the final result is 


sin 

and 
B « 

A * 


X - 2 Qn.ix-"*' = X 2 Csa.jZ" where z « x’ 


(i) 

0 - 1 )/- 


c. =1- 


B’’ 


B-A 


18I2‘' 2'« 

9,99999 99999 99990 30428 E-01 


1 143B- 55B^A 

* 18!2'" 2» 

= 2.75573 15528 52388 14908 £-06 


^ 1 , 91B^ , 13B’A 

Ciii — — TT7 T - -_ 


111 18!2* 2* 

= -2,50518 24652 10732 20541 E-08 


C,. = - 


35B' 7B^A 


13! 1812* 2’ 

= 1.60466 21504 47864 08126 E-10 


C,s = - TTT + 


B- . BA 


+ 


15! 18!2 2‘ 

* -7,35769 a3984 39792 89177 £-13 


!>■ |-1< i 

In like manner, the polynomial approximation 

. (-l)“x="*‘ 
sm X — ^ 


(2n + l)l 

can be telescoped to sin x ' x 2 dzn-iZ", 

nsv 

where 2 =x- and 

-(iy 


d. = 1- 




12 ! 2 ’= 

= 9.9999 99999 99999 783585 E-01 


C 3 = 


15B* ^ 3B^A 


_ _1 

3!"^ 18!2‘“ ■ 2“ 

-1.60666 66666 66477 96113 E-01 


I - ^4. 

d, = - r; + 


3! 1212“’ 

= -1.6606 66666 66644 563864 E-01 


C, =±- 


33B= 


21B'’A 


51 1812“ 2“ 

= 8.33333 33332 26184 74112 E-03 


^5=7:- 


5! 1212- 

= 8.3333 33332 68836 248631 E-03 


^_ aSB'^ 33B’A 

' “ 7l'‘’l8!2‘« 2» 

= -1.98412 69813 94935 49426 E-04 


d, = -i+ ™ 


7! 12!2‘ 

= -1.9841 26903 46738 823484 E-04 
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_ 1 5B- 

' ”9! 12!2^ 

« 2.7556 82887 24422 163687 £-06 


dn= - — + —--— 

' 11! 12!22 

= -2.4909 02134 88806 521004 £4)8 


and then the equations 


Eo “ Co 


at = Ci + C«bi 


<t2 “ Co+Cibi + Coba 

ag = C3+C2bi + Cjb2 + Cob3 


The remaining algorithms for sin x are all for the 
range 


for a.„ as and as. The solution is 




ao * 1 


a, = -325,523/2.283,996 


4. Fade Rational Approximations for sin x, 

|x| < - (See III.B.l.) 

6 


as = 34.911/7,613,320 

as * -479^49/11511,339,840 


Write 


bo«l 


2 Ci 2 ‘ where z=x^ and Ci ® ——— 

X iw ( 2 i + l)! 


From IV.C.2. it is known that this polynomial gives 

the desired accuracy for I x I < ^. Thus, taking 

6 

n = m = 3 in the formulas of III.B.L. 


VCi.-;bk=0 fork=4,5,6 


fork=0,1.2,3 


so that 


bx - 18,381/761.332 


b, = 1.261/4,567,992 


b. * 2,623/1,644.477,120 


ao+aiZ+a22' + aj2^ 
bo+biZ+b22* + b3Z^ 


and bo=l, 


there result the equations 


Cibs+Cabz+Cabj = 
C2b:<+C3by+C4bi — 
C3b3 + C4b2+C5bi = 


forbi, bj and bs 


Using COF 3 (Appendix A.2.a.) to find the co¬ 
efficients C,. Cl, C?, C 3 and Bi, and B 3 , diis 
rational form is now converted to continued frac¬ 
tion form 


sin x=x Co+ 


(z + Bi)+ (z+B2j+ (z+Bs 


ft 


i 


and evaluated from bottom to top as indicated in 
Appendix A.2.b. 

a = -2.6101 46506 18158 052394 EOl 

C, = 7.3922 21532 39111 816487 E03 

Ca = 7.3906 24698 46447 182801 E03 

C. = 2,4085 74343 97122 268316 E03 

B, = 1.3171 07143 04415 748255 E02 

B, = -2.5089 20716 84564 146366 EOO 

B., = 4.3867 21136 85869 652633 EOl 


and 



Solving the same system of equations as before, 
but replacing Ci by C,, the result in terms of { is 

^_ 2623 

‘ 1,235,520(1331-39901) 

23 2448b. 

■ “ 326,040 19 

b, = + 72b.: - 3024b, 

b„= 1 


Since sin x<x and iiUJL » Co + —, the right- 

X Q, 

hand side must yield a number less than one. In 

fact, when | x | < ^ , the right-hand side should be 

6 

between .955 and 1. To obtain such a number by 
adding 


to C. = -26. ... 

Q: 


ao ® 1 


.. = ({-'5)+b. 

‘■'s + O-i)'’'-"'” 


5040 


-1 . I). 


+ + b,’ 

120 V* 6/ 


could result in the loss of at least one significant 
digit. Possible remedies are to go back to the orig¬ 
inal rational form, or to inboduce a parameter { 
and an extra term (z into the computation (Ref. 5, 
p. 25) by writing 

sin X , . ^ , 

- + = 2 CizS 


where 


C, « ^ ^ —fori = 0.2,3,4.5.6 
(2i-f 1)1 


After finding at and bi for a particular value of 
the corresponding values Ci and B\ must be re¬ 
computed from COF 3. Then 


sin x«x ^Co - fz + —j . 


Values of { tried were J = .078, .0785, 2.532 and 
2.533, and the corresponding values of C© were 
about .0907, .2534, .3399 and .0145 respectively. 

^=0 is, of course, the case originally computed. 
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5« Comparison of Results 


1604 runs in double-precision using arguments from u 
lowing results: 

=0® through u=90® at 

interv'als yielded the fol- 

Method 

Maximum 
Absolute Error 

Maximum 
Relative Error 

Max. Error of N in the 
kth Significant Digit 

2. Taylor, I x I < ^ 

2.56(10-'*) at 90* 

256(10-’*) at90" 

2.56 in 16th at 90" 

3a. Cheby, | x | < g 

3.45(10-'*) at 90" 

9.64(10-'*) at 1" 

3.45 in 16th at 90' 

3b. Cbeby. | x | < g 

3.78(10-") at30* 

2.10(10-'*) atl* 

8.59 in 17th at 4“ 

4, Fade Rational, 1 x < ^ 

' ' 8 

3.57(10 ") at 30* 

7.13(10‘')at30* 

3.57 in 17th at 30" 

Cont’d Fraction, |=0 

3.57(10-") at30" 

7.13(10-") at30* 

3.57 in 17th at 30* 

" " ,078 

5.12(10 '^) at 30" 

1.02(10-'*) at30" 

5.12 in 17th at 30' 

" "1= .0785 

5.14(10 ") at30" 

1.03(10-'*) at30" 

5.14 in 17th at 30" 

" " f*2.532 

2.3304(10-") at30" 

4.661(10-") at 30' 

2.33 in 17th at 30" 

" " f=2.533 

2.3301(10 ") at 30" 

4.660(10-") at30' 

2.33 in 17th at 30" 


In these double precision runs, results for a rational function and its c'orrosponding continued fraction 
were nearly identical. 


D. TAN u 
1. Range Reduction 

Two ranges are considered, I x | and | x ]< — , 

assuming first that the argument u has already 

been reduced to the inlcn al , u I < — . 

* 2 


tan 2 

tan 3 


(t)= ^ 

(i).V2+l 


Then tan u = 


B + tan X 
1 —Btan X 


a. I X I < — 
I I - g 


b. I X I < — 
11-4 


For each u find k and x so that u = 


where k =! 0, 1, 2, 3 and 0 < x < — . 

8 


Let B =2 tan k ( — ), i.e-, 


rKi) 


+ X, 


tan 


( 1 )= 


For 


j u I < “ ) let X = u and tan u = tan x 


For 


l«l> 


let X = -u and tan u = 


tan X 
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2« Rational and Continued Fraction Forms 
Obtained from the Gaussian Continued 
Fraction (Appendix B.l) 

a. For I X I < 

I I - g 


The appropriate approximant from Appendix B.l is 



The corresponding continued fraction form ob¬ 
tained by using COF 3 of Appendix A.2.a. to find 
Cl and Bi is 


tan X = X —^ 

IbJ 

[ ao+aiZ+asZ-H-aiZ^ 1 

bo+biZ+b2Z-+b5z’+b4Z* J 

where z^x-. 

Telescoping one step would reduce the degree of 
the denominator from 4 to 3 and tlie correspond¬ 
ing continued fraction from a 4 to a 3-level one. 

A- 

Consequently we shall telescope x , or correct 
X-— , using the results of III.C.7. for an “odd” 
function. 


tan X = X 



+ 


c, Cs c, * 
(Z + B»)+ (2 + B2}+(2+B,). 


= X 




, z=x 


2 


which is evaluated via Appendix A.2.b. 


Take 



and ri 


Oi 


-1 

(2i+l) 


for i=l, 7. 


Since 


Thus, by converting the original 6 or 7 level Gaus* 
sian continued fraction to rational form, and then 
hack to another continued fraction form using 
COF 3, it has been reduced to one of three levels. 


b. For 1 X 1 < — 

I I 4 


Both tan x * x and tan x = 
B: 



S“‘>(u)= = 

2‘* 


15u^' , 45u'’ 


275u" ^ 225u‘ 189u'^ 35u^ 15u 

2 ® 2 ^ 2 ’ 2 *® 2 '* 

s, = -15/2'^ 
s, * 35/2>* 


s» = -189/2’‘ 

were tested in rational form. Since the former 

proved accurate enough, only that continued 5 ^ ^ 225/2^ 

fraction form was tried. Coefficients were found 
from Appendix A.3.a., and the continued fraction — —275/2** 

was evaluated as indicated in Appendix A.3.b. 

Six = 45/8 


3. Telescoped Rational and Continued 

Fraction Forms for I x I < 5 . (See iii.c.7.) 

4 


Si, = -15/4 
and 


From the preceding Section IV.D.2.b., we have the 
rational approximation 




— “ Sjfc+x I 


. 2 ( 7 - 


i=k 


for k = 0, 6. 
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Then 



A«+Y<,+x°(y»A<,+Y3A,+-y4A2+TtA,+y,A<) ' 
Bs + Y,X* + X*(y2 B«+YsBi + Y^Bz + Y 5B3 + YeB») 


After combining, the resulting ai and b* are 
ao = 135,135 +Yo 

= 1.3513 50000 00000 01534 86631 EOS 

a, = -17,325 + (y2+3y2+15y,+105y,+945y«) 

= -1.7336 10607 38165 56878 55239 E04 

82= 378-Y4 —IOYj —105y« 

= 3.7923 56370 39100 52361 14363 E02 

as = -1 + Y, 

= -1.0118 62505 28977 08637 24561 EOO 
bo = 135,135 

b, = -62,370+Yi + (Y2+3Y,+15Y, + 105r»+945Y,) 
= -6.2381 10607 38156 27938 77667 E04 

b2= 3150 -y,-6y.-45y5-420y, 

= 3.1549 37661 62835 53263 58151 EOS 

b, = -28 +y,+ 15y, 

= -2.8176 93975 34850 64078 19239 EOl 


and 


tan X 


' a<,+aiZ+ajZ*+a3Z* 
bo+biZ+bjz®+bjz* 


4 I ^ « 

Z = X*, X < — . 

4 


Again, coefficients and evaluation of the corre¬ 
sponding continued fraction form are obtained 
using COF 3, Appendix A.2.a.6(b. 

Co = 3.5911 01496 97721 76037 74655 E-02 
C, = -9.4381 65598 19183 41369 40110 EOO 
C2 = -1.4096 32418 00227 61516 62209 £03 
C3 = -1.5692 00421 75952 56069 73336 E02 

B, = -5.5204 04171 66464 89417 73271 EOl 
= -4.0981 70874 59656 10393 42606 EOl 
B3 = -1.5783 03284 85044 64639 80047 EOl 


4. Comparison of Results 

The following re.sults were obtained from 1604 double-preci.sion runs using arguments from u=0“ through 
11=89" at 1® intervals. 


Method 

Maximum 
Absolute Error 

Maximum 
Relative Error 

Max. Error of N in the 
kth Significant Digit 

2a. Rational & Cont'd Fraction, 

1 1 ^ 

hlSg 

tan X— xAg/Bs 

4.89(10-‘») at89® 

8.53(10-”) at89* 

4.89 in 17th at 89® 

2b. Rational & Cont'd Fraction, 




M<4 




1) tanx=xA;/B7 

4.54(10-“) at45* 

4.54(10-“) at45® 

4.54 in 16th at 45° 

2) tanx=xA8/B8 

8.71(10-“) at 45® 

8.71(10-“) at 45* 

8.71 in 19th at 45® 

3. Telescoped form of 2b (1)» 




1 1 ^ 

X < - 
1 1-4 

4.70(10-“) at45® 

4.70(10-“) at 45® 

4.70 in 16th at«® 
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E. ARCTAN u 

1. Range Reduction 


The argument 0<u<ac is reduced to one of two 
ranges: 

a. I X I <V2 ”1 

u-(V2-l) 


For 0<u<l, let x = 


l + u(\/2-l) 


and arctan u — — + arctan x 
o 


For l<u< CO. let X = 


1-u (V2-1) 
u + (\/2-l) 


Stt 

and arctan u — — arctan x. 

o 


b- I X I<tan• 


16 


ForO<u<l>sety = u, A»0. B*1. 


For l<u< 00 , set y* — , A* —. B-1 . 

u 2 

Then 


for 0<y<\/2- l»set yi = tan , «= —2- 

16 16 

and 

for V2-l<y<l, set y: = tan , «= . 

16 16 

Then 




l+yy, 

and arctan u = A+B (a+arctan x). 


2. Rational and Continued Fraction Forms 
Obtained from the Gaussian Continued 
Fraction (Appendix B.2.) 

a. I X I <V2 —1 

The appropriate approximant is 


where A,,yB|o is the quotient of two fifth-order 
poh'nomials in z, z = x-. It could be converted to a 
5-level continued fraction, but was tested only in 
its rational form because it seemed less useful than 

the methods for the range I x I <tan —^ • 
ft I I- jg 


b. i X <tan —^ . 

^ 16 

The approximants arctanx = xA^/B,; and 
arctan x=xA:/B: were tried in botli rational and 
cnntiniied fraction form. Coefficients for the con¬ 
tinued fraction forms were computed from COF 3 
and COF 4 respectively (Appendix A.2 and 3). The 
approximation arctan x - xA«/B<; was not quite 
accurate enough, hence A:/B: was telescoped (or 
A«/B,. C'orrected) as described in tlie next section to 
an approximation of the same order as Ao/B,,. 


3. Telescoped Rational and Continued 
Fraction Forms for | x | < tan ^ . 

L'sing ihc formulae in I1I.C.7. for an odd function 
and the approximants A„, Bn in Appendix B.2., and 

taking n = 6, f = tan -;^and S"''(u) = T„(u)2-'*, the 
16 

coefficients in the new R‘ are: 

a„ = 1.3513 49999 99999 99825 84406 EOS 
a, = 1.7196 24603 93687 38533 64289 EOS 
aj = 5.2490 48316 37362 32796 35437 E04 

a, = 2.2180 98888 44607 11614 67914 EOS 

b, , = 1.3513 50000 00000 00000 00000 E(K 
b, = 2.1700 74603 93685 74205 67287 EOS 
b, = 9.7799 30329 54139 12080 84660 E04 
b., = 1.0721 37452 05929 68736 47196 E04 


and 


arctan x = x 


aft+a,z + a2Z^+aflZ^' 

Lbo+biZ+bjZ*+br<z^. 


arctan x ® 



z=x^ ■xl<tan — . 
' ■“ 16 
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Coefficients of the corresponding continued frac¬ 
tion found from COF 3, Appendix A.2.a. are: 


a = 

2.0688 

56828 18530 

47509 

55450 

E-01 

c, = 

3,0086 

82092 05174 

87448 

28121 

EOO 

c, = 

-3.4976 

10177 36154 

25858 

601% 

EOO 


-1.3433 

64284 54181 

78822 

14637 

E-01 

B, = 

5.1827 

26637 17441 

95978 

52947 

EOO 

B, = 

2.6194 

66421 36919 

73145 

76466 

EOO 

B, = 

1.3197 

06666 86630 

28901 

33958 

EOO 


then 


arctan x — x 




as evaluated from Appendix A.2.b. 


4. Comparison of Results 

For arguments of the form u=tan y with y rang¬ 
ing from y»l* through y®89* at intervals of 
results of machine runs on the 1604 in double¬ 
precision are given in the table following. 



Method 

Maximum 
Absolute Error 

Maximum 
Relative Error 

Max. Error of N in the 
kth Significant Digit 

2a. 

Rational & Cont'd 
Fraction, | x | <V2 “ 1 

Arctan x=xA,n/Bio 

2.28(10-'*) at tan 45* 

4.47(10-'*) attanl* 

7.80 in 16th at tan 

2b. 

Rational & Cont'd 

Fraction, I x I < tan ^ 

16 

1) Arctan x*xA«/Bfl 

2) Arctan x = xAt/Bt 

2.42(10 '*) at tan 45* 

2.36(10-”) at tan 45* 

3.40(10-'*) attanl* 

2.75(10-'*) at tan 1' 

5.94 in 15th at tan 1* 

4.80 in 17th at tan 1* 

3. 

Telescoped form of 





2b(2),|x|<tan^ 

2.38(10") at tan 45* 

2.68(10") at tan 1' 

4.68 in 17th at tan 1* 


F. ARCSIN u 
1. Range Reductions 

No algorithm is practical for the entire range 
j Q I < 1. Arcsin u can be computed in terms of 
arctan u, or the range | u |<1 can be reduced to 

j X I < - and an algorithm for arcsin x applied in 

this range. In either case, it is necessary to take a 
square root; in the latter case, the range can be so 
adjusted that the square root operation is per¬ 
formed only a small proportion of the time. 


a. Arcsin u = arctan x where x = - _ 

quires the square root operation all the time. 

b. For 

’0<u< - , set x»u and arcsin u»arcsin X 

“ “ 2 

“<u<l, setx = 

A 

and arcsin u = — —2 arcsin x. 



re- 
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This reduction makes use of the scjuare root halt 
the time. 


c. For 


0<u< - > set x-u, A=0, B = 1 
- -2 


set x = 2ii- —1. A=, B* - 


2. Telescoped Polynomials for Arcsin x, 

lx|<i 

arcsin x _ , , ^ ^ ^ 

- 1 - ^“^6 ■^ 2 - 4-5 2 - 4 - 6-7 

1>3>5.. .39x*‘’ 

2*4-6.. .40*41 

arcsin x . 

IS the truncated series expansion of-with 

X 

error less limn 6.47(10-"‘)<2-‘". 


<u<l. 




= X- ss jso t1ial0<y<l when 0<x=< — . 


Then arcsin u = A + B arcsin x, so that a square- 
root is used .134 of the time. 


d. For 


0<u<- ,setx»ii. .A=0, B = 1 
- - 2 

1 V3 

-<u< —--sse, 

2 "2 

sctx= 2 u--l, A= •7,6= 

4 A 

— <u<--\l2+V3--965, 

2 " 2 

setx=8u'-8n- + l. A= ,B= - 


.965<u<l, 


set X® 


__iL ,A= - .B* -2 
2 2 


arcsin x 


where 


=C..+C,y + C.y^+ ... +a„y-'‘. 


Then arcsin u = A+B arcsin x and the square root 
is needed .0.35 of the time. 

All of the range reductions except a. depend upon 
an algorithm for arcsin x where I x | < - • Only one 
will be developed. 


Using the shifted Chebyshev polynomials (sec 
IU.A.2. and 3.) this imlynoinial in y of degree 20 
can be telescoped to one of degree 11. VMieii the 
transformation from the x-ariable y back to variable 
X is made, the final coefficients are; 

d„ = 9.9999 99999 99999 78634 63136 E-01 

d, = 1-6666 66666 66910 24987 33835 E-01 

d, = 7.4999 99995 41191 18650 66303 E-02 

d, = 4.4642 86051 93986 28433 58075 E-02 

d, * 3.0381 81646 51631 62166 80726 E-02 

d, = 2.2375 00912 35718 55117 84246 E-02 

d„ = 1.7312 76426 25238 66058 99121 E-02 

d, = 1.4331 24507 67095 51847 69121 E-02 

d* = 9.3428 06551 28506 27072 55181 E-03 

d. = 1.8.356 67090 64025 76498 65645 E-02 

d,„ = -1-1862 23970 78013 60943 70754 E-02 

d„ = 3.1627 12225 71360 72001 51992 E-02 
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so that 

arcsin [do + diZ+d^z^H- ... H-dnZ’*] , 

with z=X'. 

Note that one of the coefficients, diA, changed sign! 
This would suggest that perhaps the telescoping 
had proceeded too far. Such was not the case as 
results in the next section show, 

3. Test Results 

The telescoped polynomial was used in each of the 
ranges b, c, d for values of u from u ».01 through 
u= 1.00 at intervals of .01 with these results— 

Max. Absolute Erron 5.44(10*'“) atu®.5 
Max. Relative Enor: 1.04(10 **) at u=.5 

Max. Error of N in the 


3. Fade Rational 

(Diagonal of the Fade table) Ref. 18 

e* - 

P.(-x) 

where 

_ n! r " (2n-j)lxJ “I 

^ (2n)![,?u il{n-i)l J 

n=5 and n = 6 were tested. 

4. Rational and Continued Fraction Forms 
Obtained from Macon’s Even Part of 
the Gaussian Continued Fraction for e" 

(Appendix B.3.b.) 


Kth Significant Digit: S.44 in 16th at ub.S 


G, EXPONENTIAL: e- 
1. Range Reduction 

Write e"—2®e' where n and x are found as follows: 


Let y = 


log ,2 


and n » 


y ± ® the integral part 

2 


1 .. 


of y ± . 'ITie minus sign holds if u< 0 , hence y< 0 . 
2 


(Alternatively, one could compute and take the 
reciprocal when u is negative.) Let w*y—n and 
x=w Iog»2. Then e* may be computed from one of 

the algorithms which follow for |xl < ■ ^ . If 

2 

u=0, €" should be set to L 


S *^x 

e* = -where S *2+ F 

S^x 

Approximants tested for F were 

A.^ A4 Ah 

F=x 2 — , F^x* — and — 

Bs B 4 B, 

in both rational and continued fraction form. Co¬ 
efficients for the continued fraction forms were 
computed from COF 3, COF 2 and COF 2 respec¬ 
tively (Appendix A). The case F=x* A^/Bh proved 
accurate enough. 

Thus, 


6 * = 


S + x 
S — x 


2. Taylor-Maclaurin Series 

U vQ 
A=i>nl 

with 

X'" (log^)*^ 

truncation error <-< -■ ^ 1.07(10“**) 

131 2**13! 

loge2 

which is less than 2'*^ for | x | < —- 


where 


S=2+z 


2520+28z 


15,120+420z+z' 
or, in continued fraction form. 


S=2+z 


-^1 where P:=Ci(z+B 2 ) 

Q,=(z+B0(z+B.) + C. 
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and 

Ct = 28. B, = 330. 

C, = -14,580 B. = 90. 

Note that in this case the continued fraction form 
is not much shorter than the rational form. 


5. Telescoped Rational and Continued 
Fraction Forms 

Using formulae from III.C.7. for an even function, 
an unsiic'cessful attempt was made to telescope 

B, 

It would hav e rtniuced the degree of the denomi¬ 
nator from .second degree to first degree in z. 


6. Comparison of Results 


Results of 1604 double precision tests using arguments from u = —9.9 to u = 10.0 at inten als of .1 are: 


Method 

Maximum 
Absolute Error 

Maximum 
Relative Error 

Max. Error of N 
in the Kth 
Significant Digit 

3. Pad^ n*5 

3.35(10”) 

8.26(10-”) 

6,53 in 16th 

Pad^ n»6 

5.12(10”) 

1.72(10-”) 

1.33 in 19th 

4. Rational and Cont*d Fraction 




a) F = x'A,/B* 

5.76(10”) 

2.63(10-”) 

1.95 in 23rd 

b) F=x»A./B. 

5.12(10”) 

1.72( 10-”) 

1.33 in 19th 

c) Fsx’An/B, 

3.35(10 ”) 

8.26( 10”) 

6.53 in 16th 

5. Telescoped form of 4.c. 

3.33(10-*) 

3.07(10”) 

2.83 in 12th 


c 

II 

p 

at u= d:5.2 

at u=4.5 




H. LOGARITHM: Log.u = lnu 

I. Reduction of Range 


Write u=2"*m where and n may be 

zero or a positive or negative integer. 


I^t 


m-\/r/2 

X =-and compute In 

m + \/2V2 

from the algorithm. 



Then 

l„u=(n-i)l„2+l„(^) 


For 

u® 1. In u should be set to zero. 


The approxijTiations to follow are for In 
where | x |<3-2\/2. 



2. Taylor-Maclaurin Series 


, l+X ^ . . X* X* x* , 


. + 
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with error 

2x23 2-^ 10-** 


< 


•< 


23 


'7.3(10-'®) 


23 
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for 


X ,<3-2V5’. 


Achiall)’ one more term c-oiild be dropped. 


3. Telescoped Polynomial 

Let z=X"=a‘>' where a = 3—2\/2' and 0<y<l in 
the truncated Taylor series preceding. 


l+x 

c. log-- 

1—X 



in both rational and continued fraction form using 
COF 3, COF 4 and COF 4 respectively {Appen¬ 
dix A) for computing coefficients of the continued 
fraction forims. 


Coefficients for tlie continued fraction form of a, 
are: 


log 


1-i-x 

1-x 


2x 


a-y a'y- a-"y 

1 + —- -F — + ... + - 

3 5 21 


III 


Shifted Chebyshev polynomials are used and sub¬ 
stitutions made for y"’, y", >■' and y’ after which 
the coefficients are converted back to coefficients 
of z by the substitution y = 7./a-. The result h 


log 


1+x 

1-x 


2x [C. + CiZ+C,.z-+ ... +CkZ'’J ’ 


z=x= 


C., = 4,1795 91836 73469 38775 31020 E-01 

C, = -5.9412 244S9 79591 83673 46939 EOO 

C; = -3,3502 52481 31135 23355 48171 EOO 

C., = -1.2872 09952 96610 95132 66527 E-Ol 

B, = -5.1029 95.328 38691 948.33 74554 EOO 

B, = -2.5841 78755 04759 66008 33351 EOO 

B. = -1.3128 25916 56548 39157 92095 EOO 


C„ = 1.0000 000(K) (KXXX) 01720 16224 EOO 
C, = 3.3333 33333 32761 81768 85283 E-01 
C, = 2.0000 00003 09807 78908 99307 E-Ol 
C., = 1.4285 70799 46082 73472 61398 E-01 
C. = 1.1111 71831 83715 43428 06719 E-Ol 
C, = 9.0609 35658 17935 37172 14254 E-02 
C, = 8.4191 86575 86305 31375 34817 E-02 

4. Rational and Continued Fraction Forms 
Obtained from the Gaussian Continued 
Fraction (Appendix B.4.) 

.Approximations considered were 

14-X A., 

a. log- - X - 

1-x B, 

l-Fx A; 

b. log-= X- 

1-x B; 


log 


1-l-x 


1-x 


X 



5. Telescoped Rational and Continued 
Fraction Forms 

l + x A,i 

llie approximation 4.a., log --= x-— wa.s tele- 

scoped and proNided rational and a)ntinued frac¬ 
tion forms still of sufficient accuracy. The formu¬ 
las of 1U.C.7. for an odd function were used with 
n—5 and f=3 —2\/2. 


l+x 

log- 

1—x 


ao + aiz+agZ“ 

bo+b,z+b..z=+biZ' 


a,. — 

2.0789 99999 99999 84154 93231 

E04 

ai = 

-2.1545 27006 88655 98004 53920 

E04 

a. * 

4.2239 18706 18926 27409 32222 

E03 

bi. * 

1.0395 00000 00000 00000 00000 

E04 

bi - 

-1.4237 63503 44403 34724 39577 

E04 
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b, = 4.7788 37699 95350 61419 58903 £03 
b:, = -2.3041 91303 93980 93764 71785 E02 

From COF 3 {Appendix A.2.), the continued frac¬ 
tion form is 



C, = -1.8331 45841 21857 31138 11787 EOl 
C. = -2.2902 78600 16831 17096 16207 EOl 
C, = -2.3867 37346 87530 69616 24548 E-01 

B, = -15638 98338 99084 36432 27360 EOl 
B, = -3.7096 29801 61962 30733 12768 EOO 
B, = -1.3911 47833 44601 98894 76038 EOO 


6. Comparison of Results 

Letting u range from u®.! through u = 10. at intervals of .1, largest errors observed in 1604 double¬ 
precision runs are recorded below: 


Method 

Maximum 
Absolute Error 

Maximum 
Relative Error 

Max. Error of N 
in the Kth 
Significant Digit 

3. Chebyshev teles, polynomial 




a. degree 7 in z 

2.95(10-’*) 

4.25(10-’*) 

2.95 in 18th 

b, degree 6 in z 

6.24(10-’') 

9.00(10-”) 

6.24 in 17th 

4. Rational & Cont'd Fractions 




a. xA„/Br. 

6.81(10-’*) 

9.83(10-’*) 

6.81 in 18th 

b. xA:/B: 

5.11(10-’*) 

7,37(10-”>) 

5.11 in 18th 

c, xAh/Bm 

.3.83(10-»") 

5.52(10-*") 

3.83 in 20th 

5. Telescoped form of 4.a. 

7.12(10-'*) 

1.03(10-'*) 

7.12 in 16th 


u = .5,2,4,8 

u=.5, 2 

u=.5.2 
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Appendix 


Formulae for Conversion of a Quotient of Two nth 
order Pol>'nomials to Continued Fraction Form and 
for Evaluating the Resulting Continued Fraction 

1. COF 2: n = 2 

a. Conversion formulae 

aft + axZ-i-a;Z- 

LetF.^ b„+b,z+b:z- 

^ , C. C, 

(z+B,) + (z+B.) 

C)=a;5bj 

“ao“Cob<i 

aj — Bx ^Ci»bi 

Bi ® (C ibi )/tti 

C:*® (C|bo~a:oBi)/ai 

Bj ®«o/«i 


b, Evaluation 

P, ®C,(z+B,) 

Qx = (z+Bt)(z + B,)+C2 

F.=C,+ ^ 

2. COF 3: n = 3 

a. Conversion formulae 

_ ao+aiZ+aaZ*+a,z’ 
bo+b. 2 +b.z=+b,z^ 

Cl Cz Ca 

(z+B,) + (z+B,) + (z+B.) 

Co ® a^/bj 
fto ®ao~Cobo 
® Bi “Cobi 

ttg ”a2^Cpba 

Ci=af”/b3 

Bi = (biCi—a'i)/a2 

T “ byC i “ aoB I 


Cv’ = (bjCi—«,Bi“ao)/«;! 

B2 = («»C.-T)/W 
C,*(-T‘B, + «,A)/W 
B. =TAV 


b. Evaluation 

P, =C2(z+B4 

Q, «(z+B,)(z + B<)+C 3 

P, =C»Qi 

Q, = (z+B,)Qi + Pt 



3. COF 4; n = 4 

a. Conversion formulae 

ao"KaizH-a;;Z-+a;<z^+aj2* 

*■'*“ b„+b,z + b;2-+b:,z’ + b,z' 

=C„ + 

Cl C^ Ci Cj 

(2 + B,) + (z + B2) + (z + B:0 + {z + B») 

Co®ai/bi 
do ®ao^Cob<t 
ai ®ai“Cobi 
<»2 ®ai**Cob2 

Os ®a3“Cnbi 
Cl ^OTj/b* 

B, = (baCi—a2)Aa 

C 2 ® (b2Ci“Q:i'"«2Bi)/a3 

R = (aiB,+a«“biCi) 

S = (ortBj “boCi) 

W =a,C, 

B2=(tf2C2 + R)/W 

C, = (B2R+aiC2+S)/W 
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V =aw 

Bs =(C3R+a„C2+B,S)/(-V) 

C 4 = [a:oCaB3 + S(B2B3 + C3)]/( — V) 
B, = (B2S + o:oC2)/V 

b. Evaluation 

Px «C,(z+B.) 


Q. = (z+B,)(z+B4)+C, 

Pj ^CjQi 

0* = (z+B3)Q, + P, 

P3 =CiQs 

Q3 = (z+B,)Qj+P, 
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Appendix 


Gaussian Continued Fractions and Tlieir 
Approximants 

1. Tan X 

tan X 1 X® X* i* x* x* x* x* Aa 

X “ ^ ^ 7^ ^ iT^ ~ bT 


A, = l 
A.=3 
A,= 15-x2 
A3 = 105-10 x' 

A, = 945-105x2+x* 

A, = 10,395-1260x= + 21x' 

A, = 135,135 - 17.325X*+378x* - x' 

At= 2,027,025 - 270,270x2 + 6930x* -SSx* 

A,=34,459,425 - 4,729,725x2 + 135,135x‘ 
-990x‘ +x* 

2. Arctan x 

arctanx _ 1 x* 4x’ 9x* lOx* 

;; 1+ 3+ 5+ 7+ “9+“ 

Ao = l 
A,=3 

A, = 15+4x2 

A, = 105 + 55x2 

A.=945+735x2+64x‘ 

A, = 10,395 + 10,710x2 + 2079x‘ 

A, = 135,135 + 173,250x2 + 53,487x‘+2,304x‘ 

At= 2,027,025 + 3,108,105x2 + i,327,095x‘ 

+ 136,431x» 

A,=34,459,425 + 61,486,425x2 + 33,648,615x« 

+5,742,4^x«+147,456x' 

A,=654,729,075+l,332,431,100x2 + 891,080,190x* 
+2ie,602,100x' +13,852,575x‘ 

A ,0 = 13,749,310,575 + 31,426.995,600x2 

+24,861,326,490x^+7,913,505,600x* 
+865.153,575x''+ 14,745.600x‘'’ 


B„=l 

B,=3-x2 

Bt=15-6x* 

B, = 105-45x2+r* 

B, = 945 - 420x2+15x‘ 

B, = 10,3^-4725x2+210x‘ -x* 

B6=135,135-62,370x2+ 3150x<-28x« 

B, = 2,027,025 - 945,945x2+S1.975x‘ 
-eSOx'+x* 

B,=34,459.425-16,216,200x2 + 9«,945x‘ 
-13,860x‘+45x* 

25x» 36x» 49x2 64x» Six* lOOx' Aa 

11+ 13+ 15+ 17+ 19+ ~2r ~ Ba 

B„ = l 

B,=3+x2 

Bt = 15+9x2 

B, = 105+90x2+9x* 

B,=945 + lC60x2+225x‘ 

Bi = 10,395+14,175x2 + 4725x‘ +225x® 

B. = 135,135 + 218,295x2 + 99.225x^ + ll,025x« 

Bt= 2,027,025 + 3,783,780x2 + 2,182,950x‘ 
+396,900x'*+ll,025x» 

B, =34,459,425 + 72,972,900x2 +51.081,030x‘ 

+ 13,097.700x''+893,025x* 

B.=654,729,075 + 1,550,674,125x2 +1,277,025,750x‘ 
+425,675,250x*+49,116,375x’' + 893,025x'‘’ 
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B.. = 13,749,310,575 + 36,010,099,125x’ 

+34,114,830,750x‘+14,047,283,250x' 
+2,341,213,875x" + 108,05e,025x''' 


Bj=15,120+‘^’'+x‘ 

B,=332,640 + 10,080x*+42x‘ 

Bs=8,648,640+277,200x^ + 1512x‘+x‘ 


3. Exponential: e" 

a. Gaussian Continued Fraction 


_ * JL X _iL 

^ 3^ 5^ 2+ 7- 


b. Macon —even part contraction of (a) 

, , 2x x» x’' X* x^ X* 

2-X+ 6+ 10+ 14+ 18+ 22+ 26 


(2+F)+x S + x 
®'“ (2+F)^ " S-x 


where S=2+F 


, _ X* x’ X= X» x’ X* 

andF- ^ 2 + 26 

F _ 1 X* X* x’ x» 

x» 0+ 10+ 14+ 18+ 22+ 

X’ Aa 

20 ^ Ba 

A„ = l 
A, =10 
Aa = 140 + X* 

Aj=2520 + 28x* 

A,=55,440 + 756x=‘ + x‘ 

A, = 1,441,440+22,176x^+54x^ 


B.=6 

Bi=60+x* 

Bj=840+20x^ 


16x* 
9- 

B. 

Ao=2 
A.=6 

A.=30-8x* 

A,=210-110x* 

A. = 1890-1470x’ + 128x‘ 

A,=20,790 - 21.420x»+4158x‘ 

A,=270,270- 346,500x* + 106,974x‘ - 4608x" 

At=4, 054,050- 6,216,21Ox*+2,054.19Ox‘ 
-272,862x» 

A. = 68,918,850 -122,972,850x*+67,297,230x' 

- ll,484,990x‘+294,912x'' 


/ 1+X \ 

4. Logarithm: log. ) 


(S) 


2 x» 4x^ 9x^ 

1- 3- 5- 7- 


25x’ 36x» 49x* 64x’ 


11- 13- 15- 17 


B»=l 

B,=3-x* 

B,= 15-9x> 

B» = 105 - 90x*+9x‘ 

B, = 945-1050x* + 225x‘ 

B, = 10,395 -14,175x*+4725x^ - 225x* 

Ba=135,135 - 218,295x2+99,225x^ -1 l,025x« 

Bt= 2,027,025 - 3,783,780x2+2,182,950x« - 396,900x‘' 
+ ll,025x' 

B,=34,459,425 - 72,972,900x2+5i,081,030x^ 

- 13,097,700x''+893,025xs 
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i 


Constants 


log,<,e 


TT 

= 3.14159 26535 
83279 50288 

89793 

23846 26433 

log,.2 

e 

= 2.71828 18284 
71352 66249 

59045 

23536 02874 

\/35 

V2 

= 1,41421 35623 
24209 69807 

73095 

04880 16887 

V3 

= 1,7320 50807 56887 72935 27446 

v/70 

log. 10 

= 2.30258 50929 
54684 36420 

94045 

68401 79914 


log. 2 

= 0,69314 71805 
21458 17656 

59945 

30941 72321 

Off 

ten — 


= 0.43429 44819 03251 82765 11289 
18916 60508 

= 0.30102 99956 63981 19521 37388 
94724 49302 

= 1.25992 10498 94873 16476 7211 
= 5.9160 79783 09961 60425 67328 
= 8.3666 00265 34075 54797 81720 

= 0.19891 23673 79658 00691 15976 

= 0.66817 86379 19298 91999 77577 
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